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Chapter  1 


Introduction 


The  analysis  of  reinforced  concrete  structures  presents  numerous  challenges  to  the  struc¬ 
tural  analyst.  Essential  prerequisites  to  performing  effective  calculations  of  reinforced 
concrete  behavior  include 

1.  Constitutive  relations  that  accurately  model  the  damage  incurred  by  the  concrete 
during  crushing  and  cracking. 

2.  Robust  theories  and  numerical  implementations  that  can  capture  strain  softening. 

3.  An  efficient  methodology  to  model  reinforcement. 

In  previous  studies,  the  senior  author  has  investigated  various  aspects  of  devel¬ 
oping  computational  methods  for  reinforced  concrete  plate  and  shell  structures.  In  [22] 
continuum  based  plate  and  shell  theories  wre  introduced  for  finite  element  modeling  which 
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permitted  differing  constitutive  behavior  at  different  laminae  to  facilitate  representation  of 
reinforcement  through  the  thickness.  In  [11]  a  comprehensive  study  of  constitutive  mod¬ 
els  and  stress-point  integration  algorithms  was  performed.  Elastoplastic  and  viscoplastic 
theories  were  presented  and  isotropic  damage  mechanisms  included.  The  principles  upon 
which  integration  algorithms  could  be  designed  were  thoroughly  discussed  in  the  context 
of  the  convex  cutting  plane  method,  a  widely  used  procedure  in  optimization  theory  [14], 
and  application  was  made  to  the  cap  model  which  has  been  used  for  modeling  soils  and 
concrete.  In  [9]  the  state-of-the-axt  in  finite  element  modeling  techniques  for  reinforced 
concrete  plate  and  shell  structures  was  assessed  and  suggestions  for  further  research  work 
were  made.  Several  areas  delineated  as  important  in  [9]  have  been  pursued  in  this  study, 
and  elsewhere  (e.g.,  L.  R.  Herrmann  and  colleagues  axe  performing  research  on  the  bond- 
slip  problem  under  NCEL  support[5]).  A  summary  of  this  report  follows. 

In  Chapter  2  we  consider  constitutive  models  and  algorithms  for  three-dimensional 
reinforced  concrete  behavior.  We  begin  in  Section  2.1  with  a  brief  discussion  of  explicit 
and  implicit  methods  for  integrating  constitutive  models.  The  cutting  plane  algorithm 
described  in  [11]  may  be  classified  as  an  explicit  procedure,  which  entails  advantages  in 
certain  circumstances,  and  disadvantages  in  others.  The  most  proment  implicit  method, 
namely,  the  closest  point  projection  algorithm,  is  described  in  Section  2.1.1  within  the 
framework  of  a  general  elastoplastic  theory  previously  considered  in  [11].  This  algorithm 
overcomes  numerical  stability  limitations  in  the  viscoplastic  case,  at  the  price  of  somewhat 


more  involved  calculations  during  each  iteration.  Viscoplastic  regularization  has  been  pro¬ 
posed  as  a  useful  device  for  circumventing  difficulties  induced  by  strain  softening,  which 
occurs  during  concrete  crushing  and  cracking.  Consequently,  the  closest  point  projection 
algorithm  may  prove  a  useful  technique  in  the  development  of  algorithms  for  reinforced 
concrete  structures  subjected  to  severe  loading  environments. 

In  Section  2.2  we  introduce  a  class  of  rate-independent  anisotropic  elastic  damage 
models.  Previously,  we  investigated  a  simple  class  of  isotropic  damage  models  [11].  but 
the  hypothesis  of  isotropy  precludes  adequate  representation  of,  for  example,  cracking, 
which  clearly  induces  anisotropy.  The  theories  described  herein,  studied  previously  by 
Simo  [20]  in  another  context,  assume  that  the  elastic  moduli  are  damage  parameters  and 
evolve  according  to  a  damage  evolution  law.  The  failure  surface  automatically  induces 
an  anisotropic  damage  mechanism.  Likewise,  the  theory  accomodates  initially  undamaged 
anisotropic  elastic  moduli,  which  are  necessary  for  representing  even  the  linear  elastic 
response  of  reinforced  concrete  when  viewed  as  a  homogeneous  continuum. 

In  Section  2.3  we  present  algorithms  for  the  rate-independent  anisotropic  elastic 
damage  model.  Fully  implicit  algorithms  are  proposed  in  Section  2.3.1.  but  it  is  concluded 
that  these  methods  entail  calculations  at  the  stress-point  level  which  are  too  intensive 
for  practical  use  in  large-scale  finite  element  analysis.  Simpler  explicit/implicit  methods 
are  introduced  as  alternatives  in  Section  2.3.2.  The  first  of  these  algorithms  freezes  an 
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“update  direction”  at  its  initial  value  whereas  in  the  second  procedure  the  update  direction 
is  recomputed  during  each  iteration  in  an  explicit,  multi-corrector  fashion  (see  Hughes  [10]. 
Chapter  9,  for  related  ideas  in  dynamics).  The  latter  algorithm,  referred  to  as  the  updated 
explicit/implicit  algorithm,  is  viewed  as  simple  enough  for  implementatation  in  large-scale 
finite  element  programs,  while  at  the  same  time  it  “almost”  attains  full  implicitness,  a 
potential  advantage  with  respect  to  accuracy  and  stability.  As  an  example,  we  apply  the 
theory  in  Section  2.4  to  an  elliptical  failure  surface  in  pressure-deviatoric  stress  space  and 
specialize  the  updated  explicit /implicit  algorithm  to  this  case. 

A  rate-dependent  generalization  of  the  damage  model  is  presented  in  Section  2.5. 
The  classical  Perzyna  idea  of  viscous  regularization  [16]  is  invoked  in  which  the  consistency 
parameter  is  replaced  by  a  non-dimensional  switching  function,  which  turns  on  when  the 
failure  surface  function  is  positive,  divided  by  a  relaxation-time  parameter.  A  variant  on 
the  updated  explicit/implicit  algorithm  is  developed  for  the  rate-dependent  model  in  which 
a  stable  subincrementation  strategy,  as  described  in  [11],  is  employed. 

The  damage  constitutive  models  for  three-dimensional  analysis  presented  in 
Chapter  2  are  useful  for  detailed  modeling  in  regions  of  supports,  haunches,  and  transition 
zones.  They  are  also  useful  in  the  generation  and  qualification  of  plane  stress  analogs  for 
plates  and  shells  applicable  to  large-scale  structural  modeling. 
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Plane  stress  generalizations  of  the  algorithms  of  Chapter  2  are  presented  in 
Chapter  3.  The  rate-independent  and  rate-dependent  cases  are  presented  in  Sections  3.1 
and  3.2,  respectively.  Use  of  these  models  in  conjunction  with  plate  and  shell  formulations 
is  discussed  in  Section  3.3. 

In  previous  research,  attention  has  focused  on  the  cap  model  for  modeling  con¬ 
crete  (see,  e.g.,  [11]).  The  cap  model  falls  within  the  framework  of  the  general  elastoplastic 
and  viscoplastic  theories  considered  in  [11],  except  for  the  fact  that  it  is  a  particular  ex¬ 
ample  of  a  so-called  multiple  yield  surface  theory.  In  the  present  framework  of  damage 
modeling,  we  axe  likewise  concerned  with  multiple  failure  surfaces  in  analogy  with  multiple 
yield  surfaces  in  elastoplasticitv  and  viscoplasticity.  For  the  cap  model,  we  presented  previ¬ 
ously  a  concise  algorithm  for  integrating  the  constitutive  equation  and  efficiently  handling 
the  three  branches  of  the  yield  surface.  General  approaches  for  handling  multiple  yield 
surfaces  have  been  presented  in  [18,19].  Herein  we  adopt  a  more  general  point  of  view 
and  in  Chapter  4  develop  multiple  failure  surface  analogs  of  the  single  surface  models  and 
algorithms  described  in  Chapters  2  and  3.  The  approach  assumes  “m  ’  failure  surfaces. 
Consequently,  by  taking  m  =  3,  and  specifying  the  failure  surfaces  to  be  those  of  the  cap 
model,  the  approach  reduces  to  cap-like  anisotropic  damages  models.  However,  we  believe 
the  general  approach  has  considerable  potential  in  that  different  and  more  elaborate  sets  of 
failure  surfaces  will  likely  be  more  appropriate  for  modeling  reinforced  concrete  including 
cracking,  crushing,  and  bond  slip. 


The  three-dimensional  multiple  failure  surface,  rate-independent  anisotropic  elas¬ 
tic  damage  model  is  presented  in  Section  4.1.  The  updated  explicit /implicit  algorithm  for 
this  case  is  developed  in  Section  4.1.1.  Similar  developments  are  carried  out  for  the  rate- 
dependent  version  in  Sections  4.2  and  4.2.1.  Plane  stress  analogs  are  examined  in  Sections 
4.3  and  4.4.  The  rate- independent  case  is  dealt  with  in  Section  4.3  and  the  rate  dependent 
case  in  Section  4.4.  Conclusions  are  drawn  in  Chapter  5. 

The  present  work  allows  for  modeling  reinforced  concrete  by  way  of  point  con¬ 
stitutive  equations.  The  effect  of  concrete,  reinforcement,  and  their  interaction,  is,  in 
principle,  representable  by  “homogenized”,  or  distributed,  constitutive  relations.  An  ap¬ 
proach  of  this  kind  represents  an  efficient  alternative  to  those  currently  in  use  in  which 
reinforcement  and  concrete  are  modeled  separately,  and  bond-slip  is  generally  ignored. 
W  hat  is  required  to  make  an  approach  of  this  kind  a  practical  reality  is  characterization 
of  the  failure  surfaces  and  hardening  laws  for  reinforced  concrete.  All  salient  physical 
mechanisms,  e.g.,  cracking,  crushing,  bond-slip,  tension  stiffening,  shear  retention,  etc., 
are,  in  principle,  subsumable  within  such  an  approach.  On  the  other  hand,  the  practical 
and  specific  realizations  of  the  failure  surfaces  and  hardening  laws  represents  an  essential 
but  non-trivial  endeavor.  As  a  starting  point  for  such  an  endeavor  it  seems  useful  to  study 
the  failure  theory  of  fiber-reinforced  composites.  It  is  intuitevely  clear  that  fiber- reinforced 
composites  and  reinforced  concrete  have  many  features  in  common.  In  particular,  failure 
surfaces  for  fiber-reinforced  composites  have  been  studied  extensively,  e.g..  one  may  men- 
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tion  the  well-known  Tsai-Hill  and  Tsai-Wu  criteria  [121.  (An  ellipsoid  in  stress  space  is  a 
special  case  of  Tsai-Wu.  cf.  Section  2.4).  It  is  felt  that  a  worthwhile  avenue  of  approach,  tor 
developing  the  detailed  aspects  of  the  types  of  models  described  herein  begins  with  existing 
work  in  fiber- reinforced  composites.  One  of  course  would  need  to  identify  the  material  pa¬ 
rameters  of  the  model  selected.  In  principle,  this  could  be  done  by  testing  "specimens  of 
actual  fabricated  reinforced  concrete  slabs;  however,  due  to  specimen  size  and  the  number 
of  different  tests  required,  this  would  appear  impractical.  On  the  other  hand,  assuming  a 
constitutive  model  existed  for  unreinforced  concrete  that  one  had  confidence  in.  and  if  a 
sufficiently  accurate  model  of  bond-slip  behavior  existed  (see  recent  work  of  L.  Herrmann 
and  colleagues  at  U.C.  Davis[5]),  and  employing  standard  elastoplastic  modeling  of  steel 
rebars,  then  computational  tests  might  be  performed  to  determine  the  failure  surface  pa¬ 
rameters  for  the  “composite”.  These  could  subsequently  be  used  to  define  the  anisotropic 
damage  reinforced  concrete  model. 

Likewise,  failure  surfaces  used  for  unreinforced  concrete  such  as  the  well-known 
plane  stress  Kupfer  surface  [13]  —  see  Figure  1.1  —  and  its  three-dimensional  general¬ 
izations  could  be  used  as  failure  surfaces  for  anisotropic  damage  models1.  Constitutive 
equations  of  this  kind  axe  frequently  combined  with  discrete  truss-like  or  membrane-like 
equivalent  rebar  distributions  in  current  capabilities  (see.  e.g..  Cervera.  Hinton,  and  col¬ 
leagues  [2.3]  for  rep  resen  to.  tive  approaches  of  this  kind). 

lThe  failure  surtace  could  likewise  be  the  CAP  yield  surface  used  in  previous  studies 
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Figure  1.1:  Bi-axial  strength  of  concrete. 


f'c  stands  for  the  uniaxial  compressive  strength. 


An  issue  which  needs  to  be  emphasized  when  dealing  with  severely  loaded  rein¬ 
forced  concrete  is  strain  softening.  The  theoretical  and  numerical  sensitivities  engendered 
by  this  phenomenon  are  still  subjects  of  heated  discussion  in  research  circles.  Essential  in¬ 
sights  into  numerical  difficulties  (i.e.,  spurious  mesh  dependence)  brought  about  by  strain 
softening  were  presented  in  one  of  our  previous  studies  [9j.  Basically,  two  methodolo¬ 
gies  remain  in  use  to  deal  with  this  problem.  The  first  involves  specifying  key  material 
parameters  to  depend  upon  finite  element  mesh  length  scales  (e.g.,  the  plastic  modulus 
in  an  elastoplastic  theory).  This  approach  is  throughly  discussed  in  [9]  and  references 
cited  therein.  The  second  approach  employs  viscous  regularization.  Deficiencies  noted 
for  rate-independent  elastoplastic  models  used  to  represent  strain-softening  axe  avoided 
by  rate-dependent  viscoplastic  models  obtained  by  a  viscous  regularization  of  the  elasto¬ 
plastic  model.  This  approach  has  been  advocated  by  several  investigators.  Valanis  [23] 
has  provided  computational  and  theoretical  results  supporting  this  view.  One  of  the  main 
reasons  we  developed  the  rate-dependent  versions  of  the  anisotropic  damage  models  in  this 
report  was  to  provide  a  practical  means  for  dealing  with  strain  softening  behavior.  Fur¬ 
ther  research  is,  of  course,  still  necessary  on  this  topic.  Nevertheless,  within  the  spectrum 
of  methods  and  algorithms  presented  herein,  fundamental  objections  to  specific  classes 
of  models  (e.g..  rate-independent  elastoplastic  models),  may  be  overcome.  In  principle, 
the  classical  Perzyna  viscous  regularization  employed  herein  also  avoids  the  necessity  of 
employing  mesh-dependent  material  parameters  as  described  in  [9].  Likewise,  the  rate- 


9 


dependent  mechanisms  present  allow  for  more  faithful  representation  of  material  response 
to  high  rates  of  loading. 

Throughout  this  report  we  focus  our  attention  on  small-deformation  response. 
The  generalization  to  finite  deformations,  necessary  for  modeling  large  translations  and 
rotations  associated  with  reinforced  concrete  plate  and  shell  structural  response  to  severe 
blast  loadings,  may  be  simply  accomplished  by  procedures  described  in  [8].  See  also  [IS]  for 
an  update  on  approaches  of  this  kind,  and  a  description  of  their  trials  and  tribulations.  By 
way  of  these  by  now  standard  procedures,  all  that  is  developed  herein  may  be  immediately 
generalized  to  the  finite-deformation  case. 
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Chapter  2 


Constitutive  Models  and  Algorithms 
for  Reinforced  Concrete 

2.1  Explicit  and  Implicit  Methods 

In  a  previous  report  (see  [11]),  general  classes  of  inelastic  constitutive  equations  were  stud¬ 
ied  and  their  applicability  to  concrete  was  examined.  It  was  pointed  out  that  two  general 
methodologies  had  emerged  for  the  numerical  integration  of  rate- independent  inelastic 
constitutive  equations:  the  cutting  ■plane  algorithm  [14]  and  the  closest  point  projection 
algorithm  [18]. 

The  cutting  plane  algorithm  was  discussed  in  detail  in  [1 1  j.  It  avoids  the  evalua¬ 
tion  of  gradients  and  Hessians,  and  thus  is  computationally  simpler  than  the  closest  point 
projection  algorithm  and  therefore  should  be  utilized  whenever  possible.  However,  severe 
limitations  exist  when  a  viscous  regularization  is  employed:  the  cutting  plane  algorithm  is 
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explicit  and  numerical  stability  conditions  may  engender  excessively  small  subincrement  al 
time  steps  (see  [11]).  On  the  other  hand,  the  closest  point  projection  algorithm  is  im.phc.it 
and  unconditionally  stable.  Consequently,  it  offers  a  potentially  superior  alternative  to 
the  cutting  plane  algorithm  when  viscous  effects  are  present  in  the  constitutive  theory. 
However,  the  closest  point  projection  algorithm  is  more  computationally  intensive  than 
the  cutting  plane  algorithm  per  iteration,  and  thus  this  also  needs  to  be  weighed  in  any 
comparison. 

The  closest  point  projection  algorithm  will  be  illustrated  below  for  a  general 
class  of  rate-independent  elastoplastic  materials. 

2.1.1  Closest  Point  Projection  Algorithm 

The  following  notations  are  employed: 
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£pl 

/ 

q 

A 

h 

r 

c 

R 

A 

I 

TOLi 

TOL2 


Iteration  counter 
Cauchy  or  true  stress 
Total  strain 
Plastic  strain 

Elastic  strain  energy  density  function 
Yield  function 

Strain-like  internal  hardening  parameters 

Consistency  parameter 

Function  defining  the  direction  of  q 

Function  defining  the  direction  of 
Hessian  of  the  elastic  strain  energy  density  function 
Residual  vector  in  the  Newton  iteration 
Tangent  matrix  in  the  Newton  iteration 
Identity  tensor 

Tolerance  for  feasible  stress  region 
Tolerance  for  nonlinear  residual 

Partial  derivatives  are  written  in  the  following  short-hand  notation: 


«  ,  df 
a,r/  =  ^ 


dq  /  = 


d± 

dq 


d&qf  — 


02f 

dcrdq 


,  etc. 


(2.1) 

(2.2) 

(2.3) 
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The  theory  employed  has  been  presented  in  [11],  For  completeness  we  recall  it 


here.  For  further  details,  the  reader  is  urged  to  consult  [11]. 


Box  1.  A  small- deformation  rate-independent  elastoplastic  constitutive  theory. 


Constitutive  equation: 


<r  —  deyHf  —  ep^  (2.4) 

Hardening  law: 

q  =  Ah  (a,  q)  (2.5) 

Flow  rule: 

epl  =  Ar(<r,q)  (2.6) 

Loading/unloading  conditions: 


A  >  0 

(2.7) 

/(<r,q)  <  0 

(2.S) 

>-• 

£ 

II 

o 

(2.9) 
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Remarks: 


1.  Note  that  the  elastic  constitutive  law  mav  be  nonlinear  and  the  flow  rule  mav  be 


nonassociative.  In  the  associative  case 


r  =  dcrf 


(2.10) 


2.  The  loading/unloading  conditions  are  written  concisely  in  the  so-called  Kuhn-Tucker 
form  of  optimization  theory. 


3.  The  consistency  condition  yields  the  following  expression  for  A: 


A  = 


d0 rf  ■  ce 


(dof  •  cr  —  dqf  •  h) 


(2.11 ; 


The  closest  point  projection  algorithm  is  given  by  Algorithm  1  [1 7]: 


Algorithm  1.  Closest  point  projection  algorithm  for  the  small- deformation  rate-independent 
elastoplastic  theory  of  Box  1. 


Step  1.  Initialize: 


k  =  0 


p!«» 

en+l 


(2.12) 

(2.13) 
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q™i  =  q„ 

(2.14) 

«-C.  =  0 

(2.15) 

Step  2.  Compute  stress,  yield  function  and  flow  rule  and  hardening  law  residuals: 

=  Se'i  (t.+i  -  £?+,*’)  (2.1C) 

q!3.)  (2.17) 

"i+.  =  (  +  e"'  }  +  {  S'  1  (21S) 

l  -Cln  +  1  +  Qn  J  l  "n+l  J 

If  ((*  =  0  and  /<?,  <  TOLj)  or  (j/^l  <  TOL,  and  IIR^H  <  TOL2))  then 
nl  _  Dl(fc) 


Cn+1  =  *n+l 


(it) 

^n+l  —  ‘In+l 


(2.19) 

(2.20) 


return 


endif 


Step  3.  Compute  elastic  moduli  and  consistent  tangent  moduli: 

_  £Pj(^ 


Cn+1  —  d\e^  (  €n+l 


X(k]  - 

^n+l  ~ 


-71+1  J 


'(cS,)-' +  wS.a^rS,; 

-i  +  a^.aqhS,  _ 

Step  4.  Calculate  increment  of  consistency  parameter: 


fW  n  A*)  n  A* )  Ti  1 

x  x  (fc)  J  7i+ 1  -  0(Tjn+ 1  OqJn+i  An+1K.n+, 

/^n+1  _  _ 


e(k) 


(fc)  t>  (fc) 


{(fc) 

rn+  1 

i  (fc) 

1171+1 


(2.21) 

(2.22) 


(2  23) 


1G 


Step  5.  Calculate  incremental  plastic  strains  and  internal  variables: 


A*1  _  («£! ,)■'  o  w  [  ,,,  _  ^  f  rS,  y 

Aq S,  h  o'  -I  "‘P’*'  A  "‘l  C,  1J 


Step  6.  Update  plastic  strain,  hardening  and  consistency  parameters: 


pl(fc+i)  _  pl(i)  .  pl(fc) 

-n+l  —  en+l  •  <-*'«+! 


[2.25) 


a<*+»>  -  a(*)  +  a  a(fc) 

Hn+l  —  4n-t-l  +  ^4n  +  l 


(2.26) 


«!5i"  =  +  AA™, 


(2.27) 


Set  k  =  k  +  1  and  go  to  Step  2. 


Remarks: 


1.  The  above  algorithm  is  conceptually  simple:  the  rate  equations  are  integrated  via 
an  Euler  backward-difference  scheme  and  the  resulting  (nonlinear)  equations  are 
then  solved  with  Newton's  method.  A  sketch  of  the  derivation  is  provided  for  the 
interested  reader: 


Backward  difference  formulas: 

^n+1  —  ^^n+l^n+1 

(2.2S^ 

Qn  + 1  Qn  ^^n+l^n+1 

(2.29) 

=:  ^n+1 

(2.30) 
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Newton's  method: 


Linearization  is  performed  about  the  k ^  iterative  approximation  to  state  n  +  1.  The 
notation  “A”  is  used  to  denote  the  difference  between  consecutive  iterates,  e.g., 


A  \(*0  _  \(*+l) 

^An+1  —  An+1  ~  An+1 


(2.31) 


This  should  be  contrasted  with  the  use  of  u6”,  i.e., 


cX(k)  ,(*)  _A 

0An+ 1  —  An+1  Ar 


(2.32) 


We  need  to  evaluate  the  backward  difference  formulas  at  iterate  k  +  1  and  linearize 
about  iterate  k,  viz. 

e?+.‘+"  =  (2.33) 

qi+.”  =  <1„  +  aSV’kSV  (2.34) 

«S?’  +  A,*’  -  +  AAS.rS,  +  aS,  (rS,  +  Ar S.) 

=  eP‘  +  AAS,rS,+ 

<;iSi  (rSi  +  3<rrS,  •  A<rSi  +  3qrSi  •  AqSi)  (2-33) 

Likewise, 


qS,+AqS,  «  qn  +  AAS,hS,+ 


(*)  u(k) 


6Xnh  (h!2i  +  <9o-hS,  ■  Ao-15,  +  dqr^,  •  Aq^,)  (2.36) 


Linearizing  cr  results  in. 


A<r(fc)  -cW 

i-i<Tn+l  ~  Cn+1 


-  Ae 


pl(fc) 

n  +  l 


(2.37 


IS 


Substituting  this  result  into  the  previous  two.  and  combining  in  a  matrix  format 


results  in 


a,*'  ,  _ 

AqL+i 


0 


-1 


0 

-I 


aSx  (rS^  +  AA^ 


r{k) 
n+ 1 

‘n+l 


(2.3S) 


where  is  defined  in  (2.22)  and  is  defined  in  (2. IS).  In  order  to  determine 
an  expression  for  AA[,+j,  we  must  linearize  /  about  state  k ,  viz.. 


o  -/; 


n+ 1 


=  /£',+ a/: 


(*> 

+i 


/  L\ 

Solving  this  expression  for  AA„+i  yields  (2.23). 


(2.39) 


2.  The  algorithm  has  a  simple  geometric  interpretation  in  the  case  of  perfect  plasticity 
(  q  =  0  ),  linear  elastic  response  (  c  =  constant  ),  and  an  associative  flow  rule 
(h  =  dcrf).  o-n+1  is  the  projection  onto  the  yield  surface  of  the  trial  elastic  stress 
=  Se'P  ^en+1  —  taking  as  a  metric  the  tensor  c  (see  Figure  2.1).  In 
summary,  c„+i  is  the  closest  point  projection  of  onto  the  yield  surface  in  the 

energy  norm. 


3.  In  the  case  of  an  associative  flow  rule,  normality  is  enforced  with  respect  to  state 
n  +  1. 
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4.  The  algorithm  is  fully  implicit. 


2.2  Rate- independent  Anisotropic  Elastic  Damage  Model 

The  starting  point  is  the  assumption  that  reinforced  concrete  can  be  suitably  defined  by 
elastic  moduli.  We  note  that  nothing  in  the  theory  precludes  these  moduli  from  defining 
an  anisotropic  material.  For  example,  the  moduli  could  represent  the  anisotiopv  induced 
by  the  pattern  and  amount  of  reinforcement.  One  assumes  in  this  case  that  the  basic 
constitutive  relation  is  given  by 

<r  —  ce  (2.40) 

Time  differentiating  the  above  equation  one  gets 

&  =  ce  +  ce  (2.41) 

As  opposed  to  classical  plasticity,  where  a  rate  equation  defines  the  evolution 
of  the  plastic  strain,  the  proposed  damage  model  contains  an  evolution  law  for  the  elastic 
moduli  c. 

The  following  notation  is  introduced: 
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<x  Stress  tensor 

e  Strain  tensor 

c  Elastic  moduli 

a  Strain-like  internal  variable 

q  Stress-like  internal  variable  conjugate  to  a 

<p  Failure  surface 

/  Stress  dependent  contribution  to  <j> 

xp  Total  free  energy  function,  xp(e,c, a)  =  ■  ce  +  'H(a) 

H  Surface-like  energy  such  that  q  =  —H'  (a) 

A  Consistency  parameter 

<7/  Material  constant  used  in  the  definition  of  the  failure  surface 
N  Update  direction  given  by  N  =  cdcr  <P 
®  Tensor  product 

The  stress-dependent  contribution  /  to  Ahe  failure  surface  <?!>  should  not  be  con¬ 
fused  with  the  yield  surface  employed  in  the  theory  of  Box  1.  As  will  become  clear,  the 
proposed  model  contains  no  plasticity,  relying  instead  on  an  elastic  constitutive  relation 
combined  with  a  failure  surface. 

The  continuum  damage  model  is  then  given  by  the  following  equations: 
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Box  2.  A  small- deformation  anisotropic  elastic  damage  constitutive  theory. 


Constitutive  equation: 


a  =  ce 


=  /(<*■)  +  q  -oj  <  o 

A  >  0 

q)  =  0 


(2.42) 

(2.43) 

(2.44) 

(2.45) 

(2.46) 

(2.47) 
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Remarks: 


1.  Since  W  is  a  function  of  a  only,  we  employ  the  notation 

H'  =  dnH  (2.4S) 

H"  =  &aoH  (2.49) 

2.  Because  of  the  specific  form  chosen  for  the  failure  surface,  one  should  note  that 

da<P  =  daf  (2.50) 

3.  In  a  previous  approach  [11],  damage  was  defined  by  a  single  scalar  d  £  [0, 1]  that 
altered  the  elastic  moduli  in  an  isotropic  fashion.  Although  this  approach  may 
be  suitable  for  certain  applications,  the  cracking  of  concrete  intuitively  gives  rise 
to  anisotropic  damage.  The  model  under  consideration  accomodates  anisotropic 
changes  in  the  elastic  moduli  (see  equation  (2.44)). 

4.  The  total  free  energy  consists  of  a  classical  elastic  energy  and  a  surface-like  energy 
'H  that  accounts  for  damage  hardening/softening  behavior. 

5.  Unlike  classical  elasticity,  the  elastic  moduli  c  are  treated  as  progressively  degrading 
internal  parameters;  these  are  determined  through  their  initial  values  —  possibly 
ar.’sotropic  to  account  for  reinforcement  —  and  an  evolution  law  which  redefines 
them  in  specified  directions  according  to  the  failure  ciitcrion.  An  illustration  of  this 
behavior  is  shown  in  Figure  2.2. 
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Figure  2.2:  Stress-strain  relation  for  a  one- dimensional  elastic  damage  model.  The  damage 
mechanism  causes  unloading  at  reduced  values  of  the  elastic  modulus.  Unloading  occurs 
towards  the  origin. 
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6.  The  equations  in  Box  2  can  be  obtained  by  employing  the  principle  of  maximum 
dissipation  (see  [1.6. IS]).  We  sketch  the  derivation  below. 


Time  differentiating  the  free  energy  leads  to 


rif  =  €  •  ce  4-  -e  •  ce  +  Ti'a 

9 


(2.51) 


By  definition,  the  dissipation  function  T>  gives  the  rate  of  change  in  free  energy  due 
to  changes  in  the  inelastic  internal  variables.  From  (2.51),  we  get 


T>  =  —  -e  ■  ce  —  Ti! a 


(2.52) 


Hence,  the  rate  of  change  in  the  free  energy  is  given  by 


y>  =  ce  —  V 


(2.53) 


The  second  law  of  thermodynamics  requires  that,  for  all  e, 


— if  +  (T  ■  e  >  0 


(2.54) 


Satisfaction  of  (2.54)  can  only  be  obtained  if  the  following  two  conditions  hold 


cr  =  ce 


(2.55) 


V  >  0 


(2.56) 


Hence,  the  second  law  of  thermodynamics  combined  with  the  form  of  the  free  energy 
assumed  leads  to  two  conclusions: 


(a)  The  stress-strain  relation  takes  on  the  usual  elastic  form,  namelv  cr  —  ce. 
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(b)  The  dissipation  function  P  is  non-negative. 


In  order  to  state  the  principle  of  maximum  dissipation  we  introduce  the  set  of  ad¬ 
missible  states  £, 

£  =  {((T,q)  \<p{(T,q)  <  0}  (2.57) 

The  principle  of  maximum  dissipation  then  states  that  the  dissipation  function  P  is 
maximum  subject  to  the  constraint 


{tr,q)  e  £ 


(2.58) 


Equivalentv,  this  can  be  phrased  in  terms  of  a  minimization  problem  for  — P,  namely 

(-P) 


mm 

(<r,q)  €  £ 


(2.59) 


If  it  were  not  for  the  constraint  <P  <  0.  the  optimality  conditions  could  be  simply 
obtained  by  differentiating  — P  with  respect  to  <r  and  q.  This  approach  can  still  be 
utilized  if  one  employs  the  method  of  Lagrange  multipliers,  which  transforms  the  con¬ 
strained  minimization  problem  above  into  an  unconstrained  problem  by  appending 
to  — P  the  constraint  times  a  Lagrange  multiplier  A.  This  gives  rise  to  a  Lagrangian 
£,  defined  by 


£  =  -P  +  Ao  (<r.  q) 


(2.60) 


The  optimality  conditions  are  then  given  by 

Oa-C  —  0 
0qC  =  0 


(2.61) 

(2.62) 
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along  with  the  Kuhn-Tucker  conditions 


A  >  0  (2.63) 

0(tr,q)<  0  (2.64) 

\<t>  (<r,  q)  =  0  (2.65) 

In  order  to  perform  the  differentiations  in  (2.61)  and  (2.62),  we  rewrite  the  dissipation 
function  T>  as  an  explicit  function  of  <r  and  q ,  obtaining 

C  =  •  c_1cc_1<r  -  qa  +  A  [/  (er)  +  q  —  07]  (2.66) 

Differentiation  with  respect  to  <r  and  q  is  now  straightforward  leading  to 

ce  =  -Xcdaf  (2.67) 

d  =  A  (2.68) 

Note  that  the  consistency  parameter  A  can  be  interpreted  as  a  Lagrange  multiplier 
that  forces  <r  and  q  to  stay  within  the  allowable  region. 

7.  In  view  of  (2.68),  the  hardening  law  (2.43)  can  be  rewritten  as 

q  =  -*dlaH  (2.69) 

8.  During  amage  loading,  o  =  0.  From  (2.45)  one  gets 

<Z>  =  d<rf  ■  &  +  q  (2.70) 


2S 


Time  differentiating  the  constitutive  equation  (2.42)  and  substituting  into  (2.70) 
leads  to 


o  =  dcrf  ■  {ck  +  ce)  +  q 


(2.71) 


Substituting  the  hardening  and  damage  evolution  laws,  equations  (2.43)  and  (2.44), 
into  (2.71)  and  setting  <p  to  zero  results  in  an  explicit  expression  for  the  consistency 
parameter, 


dcrf  ■  ck 

dcrf  ■  cdcrf  +  'H" 


(2.72) 


9.  The  derivation  of  the  expression  for  c,  (2.44),  deserves  special  attention.  Assuming 
that  c  is  a  rank-one  tensor,  the  symmetry  conditions 


Cijkl  —  cklij  —  Cji  kl  —  C\]lk 


(2.73) 


and  (2.67)  imply  the  damage  evolution  law: 

c  -  )  cd<r/  $  cdg-f 

dcrf  •  ce 


(2.74) 


The  rank-one  assumption  is  the  key  to  obtaining  this  simple  expression.  Note  that 
(2.67)  and  (2.73)  can  be  satisfied  by  infinitely  many  other  definitions  of  c.  This 
may  be  a  worthwhile  avenue  of  future  research  in  that  cracking  and  crushing  may  be 
amenable  to  more  realistic  treatment  by  appropriate  generalizations  of  the  definition 
of  c.  The  subject  of  quasi-Newton  updates  may  be  relevant  in  this  regard  (see.  e.g.. 
Luenberger[l4j)  in  that  (  2.67)  has  the  form  of  the  so-called  quasi-Newton  equation.  It 
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should  be  possible  using  the  methodology  of  quasi-Newton  updates  to  satisfy  (2.G7), 


(2.73),  and  other  attributes  deemed  useful  in  the  descrpti.  :i  of  cracking,  crushing, 
bond-slip.  etc. 

2.3  Algorithms  for  the  Rate-independent  Anisotropic 
Elastic  Damage  Model 

2.3.1  Fully  Implicit  Return  Mapping  Algorithm 

A  first  attempt  to  integrate  the  above  equations  might  employ  a  similar  approach  to 
that  used  in  the  closest-point  projection  algorithm  for  plasticity.  This  would  lead  to  the 
following  set  of  nonlinear  equations: 

c  \  fn+l  0  &GT  fn- 4*1  /  r.  —  «-  \ 

Cn+i  —  Cn  CM  —  -  (2.10) 

Oct J  ■  cn+ien+i 

crn+1  =  cn+,en+i  (2.76) 

qn+l  =  -H(an  +  6\)  (2.77) 

In  an  elastic  step,  one  has  6\  =  0  and  the  solution  of  the  above  is  obviously  trivial. 
However,  if  damage  evolution  occurs,  to  the  above  equations  one  adds 

<i>((rn+x,qn+])  =  0  (2.7S) 

Equations  (2.75)-(2.7S)  amount  to  a  fully  implicit  nonlinear  system.  One  has 
to  solve  yimultane.ounly  for  cn+i.  <rn+1,  qn+l.  and  6 A.  Solving  such  a  large  nonlinear  sys- 


tem  of  equations  at  each  integration  point  at  each  nonlinear  iteration  of  a  finite  element 


analysis  appears  impractical  and  alternative  integration  schemes  are  thus  called  for.  Some 
simplifications  are  described  below. 


2.3.2  Explicit/Implicit  Methods 


We  consider  the  explicit  integration  of  the  elastic  moduli  while  mantaining  an  implicit 


integration  for  the  other  variables.  This  leads  to  the  following  system  of  equations 


-trial  „  ^ 

trn+ 1  —  cnen+l 

(2.79) 

C  drr  f trial 

_ _ n'-nW/  + 1  W  CnGVr/n  +  1 

Cn+ 1  —  cn  vA  .  •  i 

Jo-fnl ld  •  c„en+i 

(2.80) 

&n+ 1  =  Cnt-l'.+l 

(2.S1) 

9n+ 1  —  —W  {otn  +  SX) 

(2.S2) 

4>  (<Tn+  li  <7n+l)  —  0 

(2. S3) 

The  above  system  can  be  solved  by  the  following  algorithm: 
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Algorithm  2.  Explicit/implicit  algorithm  for  the  small-deformation  elastic  damage  con¬ 
stitutive  theory  of  Box  2. 


Step  1.  Initialize: 

a(0)  =  0  (2.84) 

o-S,  =  ^  =  C„en+1  (2.S5) 

N™  =  c„  ■  $,/£!.  (2.86) 

k  =  0  (2.87) 

Step  2.  Update  stress  and  hardening  parameters: 

<2-88) 

<*!2,  =  a,  +  <SAW  (2.S9) 

9^1  =  (<$,)  (2-90) 


-(*)  _  r(fc)  Ifc) 

^n+1  /n+ 1  9n-t-l  ® ] 

If  ((it  =  0  and  <p[kU  <  TOL  )  or  (|^+il  <  TOL  ))  then 


(2.91) 

(2.92) 

(2.93) 
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an+l  —  an+l 


(2.94) 


return 

endif 

Step  4.  Compute  <$A-increment 

WS,  =  -  [d<rCx  •  c„N +  H"  (a™,)]  (2.05) 

6\{k+l)  =  6\(k)  -  1  (2.96) 

Set  k  =  k  +  1  and  go  to  Step  2. 

Remarks: 

1.  The  trial  stress  =  o^+i  is  similar  to  that  introduced  in  plasticity:  the  inelastic 

parameters  are  frozen  at  their  previous  values  and  the  stress  is  updated  with  the 
current  strain  value. 

2.  Equation  (2.88)  is  derived  as  follows: 

^n+1  =_  Cn+l^n+1 

—  Cn^n+l  (Cn+1  Cn)^n+1 

=  ^n+f1  +  (Cn+1  ~  Cn)fn+1  (2.97) 
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The  second  term  is  replaced  by  a  discrete  version  of  (2.67).  namely 

(c»+i  -  cn)  e„+i  =  -SX Cnda-fn+i  (2.98) 

3.  The  above  system  of  equations  contains  a  single  scalar  unknown  S A.  a  considerable 
simplification  compared  with  the  fully  implicit  system  described  previously.  This 
single  scalar  unknown  is  the  solution  of  (2. S3),  which  can  be  written  in  the  form 

0  =  4>(6\)  =  f  (<7^  -  -  H'(an  +  6X)  -  Of  (2.99) 

Step  4  of  the  algorithm  amounts  to  a  Newton  iteration  method  for  solving  this 
nonlinear  scalar  equation. 

4.  Note  that  the  update  direction  is  held  fixed  at  its  initial  value  N^.  In  an  attempt 
to  improve  upon  the  accuracy  of  the  explicit /implicit  solution  described  above,  we 
consider  a  variant  below  which,  at  each  nonlinear  iteration,  recomputes  the  update 
direction  N  =  c dcr f  ■ 
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Algorithm  3.  Updated  explicit/implicit  algorithm  for  the  small-deformation  anisotropic 


elastic  damage  constitutive  theory  of  Box  2. 


Step  1.  Initialize: 


<5  A(0)  =  0 


(2.100) 


__(°)  _  -trial  _ 

°^n+ 1  —  ^n+1  cnen+l 


(2.101) 


_  r  n  f(°) 

^  '  n+ 1  —  cnO(TJn+ 1 


(2.102) 


c(0)  -c 

cn+l  ~  cn 


(2.103) 


k  =  0 


(2.104) 


Step  2.  Update  stress  and  hardening  parameters: 


W  _  trial  c\(*)T\T(fc) 

n+l  —  "n+l  0A  iNn+l 


(2.105) 


a^i  =a»  +  #> 


(2.106) 


„(*)  _  v  f-W  A 
9n+ 1  't  J 


(2.107) 


Step  3.  Check  for  failure  and  convergence: 


A*)  _  Ah)  .  (*) 

<>>n+ 1  —  /n+l  +  9n+l  —  °7 


(2.108) 


If  ((fc  =  0  and  4>nli  <  TOL  )  or  (K+il  <  TOL  ))  then 


(*) 

an+  1  —  &  n+  1 


(2.109) 
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(2.110) 


Cn+1  —  Cn+, 

„  _  Jk) 

^n+1  ^-n+1 


return 


endif 


Step  4.  Compute  <SA-increment: 

wi+,  =  -  [So-zS,  •  n!2.  +  n"  K‘+.)] 

<A'‘+"  =  SX  '*>  - 


Step  5.  Update  return  direction  and  moduli: 


n!S-,'>  =  cS,a^/«‘> 

c(fc+!)  _  c  S\tk+l)^'n+l  *  ®  } 


sn+ 


.  N(/c+1) 

1  ^n+l 


Set  k  =  k  +  1  and  go  to  Step  2. 


(2.111) 


(2.112) 

(2.113) 


(2.114) 

(2.115) 


Remarks: 


1.  Despite  the  explicit  update  in  the  return  direction  N,  the  algorithm  is  implicit  with 
respect  to  the  solution  of  the  (nonlinear)  equation 


0  (o’ n+u  Qn+l  )  —  0 


(2.1 1G) 
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Analogous  to  (2.99),  we  write 


0  =  0  (8X)  =  f  -  6X -H'{an  +  6\)-<Tf  (2.117) 

Equation  (2.112)  follows  from  (2.117)  by  employing  (2.105)  and  (2.106). 

2.  The  update  for  the  return  direction,  equation  (2.114),  is  an  attempt  to  approximate 
the  implicit  algorithm  sketched  in  Section  2.3.1. 

3.  An  important  distinction  between  the  present  algorithm  and  the  implicit  algorithm 
of  Section  2.3.1  is  that  here  equation  (2.75)  is  not  necessarily  satisfied  by  cn+i. 

4.  The  present  algorithm  represents  an  attractive  balance  between  computational  effort 
and  implicitness.  It  is  employed  as  a  basis  for  subsequent  developments  within  this 
report. 

2.4  An  Example:  Elliptical  Failure  Surface  in  Pressure- 
Deviatoric  Space 

As  an  example  of  a  typical  application,  we  consider  a  failure  surface  given  by 


f(<r)  =  j^(p-Po)2  +  ||dev<r]|2  -  R 2 

(2.1  IS) 

p  =  — tr<r/3 

(2.119) 

Qer  =  dev<r  =  <x  +  pi 

(2.120) 

where  R,  5.  and  po  are  material  parameters. 
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Remarks: 


1.  Q  is  a  projection  operator  that  extracts  the  deviatoric  part  of  a  symmetric  tensor. 
That  it  is  indeed  a  projection  can  be  seen  from  the  fact  that 

QV  =  Q<r  (2.121) 

Along  with  Q  we  can  introduce  another  projection  operator  P, 

Pa  =  tr  a\  (2.122) 

O 

where  I  stands  for  the  identity  matrix.  Note  that 

Pa  +  Qa  —  a  (2.123) 

and  that 

PQ  =  QP  =  0  (2.124) 

2.  The  above  failure  surface  represents  an  ellipse  in  p  —  deva  space,  with  R  being  the 
semi-axis  of  the  ellipse  in  the  deviatoric  plane  and  5  the  semi-axis  in  the  pressure 
plane,  as  shown  in  Figure  2.3. 

3.  As  mentioned  in  [11],  if  the  cutting  plane  algorithm  is  employed,  the  failure  surface 
has  a  preferred  representation  which,  in  the  current  case  is 

/(<>■)=  y'f](p-po)!  +  l|devof (2.125) 


3S 


4.  po  is  introduced  in  order  to  control  the  tensile  resistance  of  the  concrete  model. 

5.  The  failure  surface  above  is  isotropic.  A  more  general  surface  of  the  Tsai-\\u  type 
(see.  e.g.,  Jones  [12])  is  given  by 

/(*■)=  |)(p  -  P»)!  +  Iqt  •  MQ<t  -  A2  (2.126) 

where  the  tensor  M  could  account  for  anisotropic  elastic  behavior  as  that  emanating 
from  reinforcement.  For  simplicity  in  the  exposition,  we  proceed  with  the  surface 
described  by  (2.118).  The  generalization  to  (2.126)  is  straightforward. 

Damage  will  be  described  according  to  the  relation 

7i  (a)  =  —Ha  (2.127) 

where  H  is  assumed  to  be  a  positive  constant.  The  negative  sign  is  indicative  of  the 
presence  of  softening  behavior  attributable  to  crushing  and/or  cracking  of  concrete. 


For  a  material  characterized  by  (2.118)  ,  one  can  utilize  Algorithm  3  in  which 


d<r/  =--~(p-Po)I  +  2Qo- 

(2.12S) 

II 

1 

bq 

(2.129) 

o 

II 

£ 

(2.130) 

For  the  sake  of  completeness,  Algorithm  3  is  specialized  to  this  case. 
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Figure  2.3:  Elliptical  damage  surface  in  the  pressure- deviatoric  plane. 
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Algorithm  4.  Updated  explicit /imp licit  algorithm  for  the  small- deformation  anisotropic 
elastic  damage  constitutive  theory  of  Box  2  with  assumed  elliptical  failure  surface  in  pressure- 
deviator  space  and  constant  hardening  modulus. 


Step  1.  Initialize: 


a(0)  =  o 


_(°)  _  -trial  _  _  , 

an+ 1  —  "n+l  ~  Cn£n+1 


JO)  _  tr_(0)  /o 

Pn+ 1  —  ~tr<Tn+l  /  3 

,  (o)  (0)  ,  (0)  T 

devo-^1  =  <7^1  +PnU1 


_  Or 

^n+1"“CnL  3  52 
(0) 

C„+l  =  Cn 


1  R 


-  (pSi  “  Pd)  I  +  devo- 


(0) 

n+1 


k  —  0 


Step  2.  Update  stress  and  hardening  parameter: 

*i+i  =  -  *A(fc)N£, 

Pnl  1  =  -toVn+l/Z 
devcr^,  =  cr{nkl,  +  p^,I 


(2.131) 

(2.132) 

(2.133) 

(2.134) 

(2.135) 

(2.136) 

(2.137) 


(2.138) 

(2.139) 

(2.140) 

(2.141) 
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Step  3.  Check  for  failure  and  convergence: 


«£.  =  ^(pS,  -  P«)'  +  ||<iev«rS,||2  -R‘  +  H  -a,  12.1421 

If  ((Jfc  =  0  and  <p(k\x  <  TOL  )  or  (1^,1  <  TOL  ))  then 


(Jfc) 

n  =  <+1 

(2.143) 

r  -  r{k) 

cn+l  —  Cn+1 

(2.144) 

(Jfc) 

«n+l  =  an+1 

(2.145) 

return 

endif 


Step  4.  Compute  (SA-increment: 


i  =  -2 


’  1  R2 

(  (Jt)  \ 

3  52 

\P»+1  -  PO) 

iNn+l 


«5A(4+1)  =  6\{k)  -  4>(klx/D<t>(k) 


n+l 


Step  5.  Update  moduli  and  return  direction: 


Ni++i1)  =  2cSi  “  ~  (p!,+j  -  Po)  I  +  devo-JS, 

c(^)_c  a(,+i)CLi^ 

C»+*  -Cn-^^  ,k  +  x) 

^n+1  *  *^n+l 

Set  k  =  k  4-  1  and  go  to  Step  2. 


(2.146) 


(2.147) 


(2.14S) 

(2.149) 
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2.5  Rate-dependent  Damage  Model 


The  rate-dependent  damage  model  can  be  obtained  from  the  rate-independent  model  by 
means  of  a  viscous  regularization  of  the  Perzvna  type[l  1.16] .  In  this  generalization.  A  is 
replaced  by  (\  (<p))/r,  where  \  is  a  non-dimensional  function  of  o.  (.)  denotes  the  Macaulay 
bracket,  viz. 


(x) 


x  If  x  >  0 
0  otherwise 


(2.150) 


and  r  is  the  relaxation  time.  Also,  the  damage  loading/unloading  conditions  are  dispensed 
with.  As  will  be  shown  below,  this  model  has  the  advantage  of  being  trivially  implemented 
in  a  finite  element  context. 


The  theory  is  summarized  in  Box  3. 


Box  3.  A  small- deformation  anisotropic  rate -dependent  damage  constitutive  theory. 


Constitutive  equation; 


<r  =  ce 


Hardening  law: 

g  - - H 

T 


(2.151) 


(2.152) 
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Damage  evolution  law: 


(\(o))  c dgrf  3  cdg-f 
t  daf  ■  ce 


(2.133) 


Remarks: 


1.  Just  as  in  viscoplasticitv,  the  function  x  may  be  taken  to  be 
X(^)  =  (sgn  <t>)(<p/y)N 
where  y  and  N  are  positive  constants,  [11]. 


(2.154; 


2.5.1  Algorithm  for  the  Rate-dependent  Damage  Model 


Integration  of  the  equations  in  Box  3  can  be  done  according  to  the  following  algorithm: 


Algorithm  5.  Updated  explicit /implicit  algorithm  for  the  small- deformation  anisotropic 
elastic  damage  constitutive  theory  of  Box  2 ,  or  rate -dependent  damage  constitutive  theory 
of  Box  3. 


Step  1.  Initialize: 


J  true,  rate-dependent  damage 
(  false,  elastic  damage 


(Box  3) 
(Box  2) 
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lquit  =  false 

-A  d.  +  1  —  ln 

(2.155) 

<5A,0)  =  0 

(2.156) 

_  t 

ln+l  —  ln 

(2.157) 

(o)  trial  _  , 

^n+l  —  <Tn+ 1  —  Cid+1 

(2.158) 

N<0)  —  r  dtr  f{0) 

CnC/CT /n-fl 

(2.159) 

(0) 

Cn+l  =  Cn 

(2.160) 

O 

II 

-Si 

(2.161) 

o-S,  =  -r'jf1  -  a«‘»NS, 
a!5l  =  a„  +  i\,t} 


(2.162) 

(2.163) 

(2.164) 


Step  3.  Check  for  failure  and  convergence: 

jl(*)  _  Ak)  ,  Jk) 

^n+l  _  Jn+ 1  +  Qn+\  ~  °  1 

If  ((Jk  =  0  and  <  TOL  )  or 
(  lvisc  =  false  and  <  TOL  ))  then  lquit  =  true 

If  (  lquit  =  true  )  then 

_  -  rrik) 

<Tn+i  —  <Tn+I 


(2.165) 


(2.166) 
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„  _  Jk) 

Cn-fl  ^n+1 


(2.167) 


^ 

an+l  —  an+l 


(2.168) 


return 


endif 


Step  4.  Compute  6A-increment: 


Dtf+i"  =  -  [At/'?,  ■  Ni*l,  +  H"  (ttS,)] 


(2.169) 


If  (  lvisc  =  false  )  then 


=  6\W  -  4>nll! D<t>nl\ 


(2.170) 


~{k)  ‘ 

in+' =  "xS’.w!3. 

A«S,  =  ti+,  1  -  exp 
,(*+l)  _  Ak)  A  Ak) 

cn  +  l  —  cn+l  T  AAln+l 


(2.171) 


(2.172) 


(2.173) 


^  (*n+i  -  in+i  >  0)  then 


At{k]  -  At{k)  -  (t{k+l)  f  ) 

^n+X  —  L\Zn+ 1  —  (tn+i  —  tn+ 1  j 


(2.174) 


lquit  =  true 


endif 


«lt+»  =  «<*>  +  A^.v^./r 


(2.175) 


endif 
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Step  5.  Update  return  direction  and  moduli: 


_  1*1  a  ft*) 

x>ln+l  —  Cn+\°(T  Jn-f\ 


(fc) 


c(k) 


(k+l) 
n+ 1 


Cn+1  -Cn  ^  —  rpiTTF 

€n+\  iNn+l 

Set  k  =  k  4-  1  and  go  to  Step  2. 


(2.176) 

(2.177) 


Remarks: 

1.  As  discussed  in  [11],  the  approach  taken  herein  obviates  the  issue  of  stability  by 
automatically  defining  a  priori  stable  sub  incremental  time  steps  within  the  context 
of  the  updated  implicit/explicit  algorithm.  The  reader  is  encouraged  to  note  the 
similarities  between  the  subincrementation  strategy  described  in  Algorithm  5  and 
the  one  in  [11],  as  applied  to  the  cutting  plane  algorithm. 
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Chapter  3 


Plane  Stress  Models  —  Applications 
to  Plates  and  Shells 


The  application  of  the  damage  models  described  previously  leads  to  a  trivial  implementa¬ 
tion  in  three-dimensional  geometries.  Evidently,  one  might  attempt  to  model  thick  shells 
and/or  plates  with  three-dimensional  elements;  however,  as  the  shell  gets  thinner,  the 
cost  of  three-dimensional  analysis,  and  numerical  difficulties  engendered  by  thin,  three- 
dimensional  elements,  suggest  the  use  of  shell  structural  elements. 

3.1  Plane  Stress  Algorithm  for  the  Rate-independent 
Anisotropic  Elastic  Damage  Model 

A  common  feature  of  the  structural  plate  and  shell  elements  is  the  zero  normal  stress  (ZNS; 
constraint  which,  in  the  two-dimensional  case,  is  equivalent  to  a  plane  stress  formulation. 


4S 


A  general  approach  for  the  analysis  of  plane  stress  problems  is  based  upon  the 
introduction  of  e33  as  an  independent  variable,  where  the  third  direction  is  normal  to  the 
plane  of  interest.  Solution  of  the  resulting  nonlinear  system  of  equations  can  be  obtained 
by  introducing  the  ZNS  constraint  as  an  additional  condition  to  be  satisfied: 

<r33  =  0  (3.1) 

One  can  easily  modify  Algorithm  3  by  noting  that  the  ZNS  constraint  can  be 

written  as 

(c^S1)  •  e33  -  a(fc+1> N&  •  e33  =  0  (3.2) 

where  633  stands  for  the  unit  tensor  in  the  33-direction,  namely 
'0  0  0' 

633=  0  0  0  (3.3) 

0  0  1 

Hence,  the  implicit  computation  of  is  done  in  conjunction  with  the  com¬ 

putation  of  e33+1*.  A  complete  plane  stress  algorithm  is  then  given  as  follows: 

Algorithm  6.  Updated  explicit/implicit  plane  stress  algorithm  for  the  small- deformation 
anisotropic  elastic  damage  constitutive  theory  of  Box  2. 

Step  1.  Initialize: 

<5A(0)  =  0  (3.4) 
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€  =  *n+l 


(3.5) 


6(0)  -  — 
e33  — 


cn.33116ll  +  cn.3322e22  +  2cn.33i2ei2  +  2cn,3323623  +  2cn,33i3fi3 


cn,3333 


_(o)  _  -trial  _ 


^n+l  —  ^n+l  —  Cn£n+1 


N<°>  =  cnd<rf, 


'n+l 

c(°)  _c 
cn+l  —  Ln 

k  =  0 


(0) 

n+l 


(3.6) 

(3.7) 

(3.8) 

(3.9) 
(3.10) 


Step  2.  Update  stress  and  hardening  parameter: 


o-Si  =  C „«£,  -  «<*> NiS, 

Q<«l  =  an  +  «'‘) 

KM 


Step  3.  Check  for  failure  and  convergence: 
aW  _  fW  ,  lk) 

<Pn+ 1  =  /n+l  +  9n+l  “  °  ! 

If  (( k  =  0  and  <  TOLi)  or 

(|*!2il  <  TOLj  and  \a{n%x\  <  T0L2))  then 
(fc) 

<Tn+l  =  <T„+ 1 
(k) 

*-n+ 1  ^"n+1 

_ 

*-*n+l  *-*n+l 

return 


(3.11) 

(3.12) 

(3.13) 


(3.14) 


(3.15) 

(3.16) 

(3.17) 
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endif 


Remarks: 

1.  Except  in  the  33-component,  which  is  updated  to  satisfy  the  ZNS  constraint,  e^1* 
is  identically  equal  to  e  =  €n+1, 

2.  Even  though  the  traditional  plane  stress  analysis  does  not  contain  the  transverse 
shear  stresses  and  strains,  they  were  included  in  the  analysis  here  so  the  algorithm 
can  be  applied  to  thick  plates  and  shells,  where  these  quantities  are  important. 
Likewise,  the  preceding  algorithm  may  be  reduced  to  the  traditional  case  by  omitting 
the  transverse  stresses  and  strains. 
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3.  The  algorithm  allows  for  different  tolerances  in  the  satisfaction  of  each  of  the  con¬ 
straints.  namely  TOLi  for  the  failure  surface,  and  TOL?  for  the  ZNS  constraint. 


4.  The  2x2  system  of  equations  one  solves  in  Step  3  reflects  the  satisfaction  of  the  two 
scalar  constraints. 

d,  (c^y,"  -  -W'  (q„  +  = 

kI(«i‘+,,,4,yi,j3)  =  *l‘+”  =  o  <3-21) 

(c,4*++,11)  ■  *33  -  «l‘+,>Ni2,  ■  e33  = 

*,  -  *!•*•»  -  0  (3-22) 


Standard  apphcation  of  Newton's  method  to  the  above  two  equations  leads  to 


ds\hi 
ds*h  2 


1.33^1 
^n+ 1,33^2 


a<fc+1>  -  a<fc) 

Jk+ 1)  Jk) 
cn+l,33  en+l,33 


(3.23) 


Differentiating  (3.21)  and  (3.22)  we  get 


dsxhi  =  -  {daf^x  •  NS,  +  H"  (<*&)) 

dln+u33hx  =  [cnd<r/i+ 1]33 


(3.24) 

(3.25) 


d$\h2 


-N 


r(k) 


n+1,33 


(3.26) 


^tn+1,33  ^2 


Cn,3333 


(3.27) 


Combining  (3.23)  to  (3.27)  leads  to  the  final  expression  in  Step  4.  equation  (3. IS). 


3.2  Plane  Stress  Algorithm  for  the  Rate- dependent 
Damage  Model 


In  this  Section  we  present  an  algorithm  which  incorporates  the  plane  stress  constraint  in  the 
rate-dependent  damage  constitutive  theory.  Our  approach  combines  features  of  Algorithms 
5  and  6.  In  the  absence  of  rate-dependent  effects,  the  algorithm  reduces  identically  to  the 
plane  stress  elastic  damage  case,  namely  Algorithm  6.  The  procedure  is  given  as  follows: 


Algorithm  7.  Updated  explicit/implicit  plane  stress  algorithm  for  the  small- deformation 
anisotropic  elastic  damage  constitutive  theory  of  Box  2 .  or  rate- dependent  damage  consti¬ 
tutive  theory  of  Box  3. 


Step  1.  Initialize: 


lvisc  = 


true,  rate-dependent  damage  (Box  3) 
false,  elastic  damage  (Box  2) 


Iquit  =  false 


'Afn+l  —  1  —  t 


n+1  Ln 


d°l  —  i 

cn+l  —  lr 


(3.2S) 

(3.29) 


6  A(0)  =  0 


(3.30) 


€  =  fn+l 


(3.31) 
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(0)  _  Cn.3311^11  +  C«,3322e22  +  2cn,33i2ei2  +  2cn,3323£23  +  2cn,33i3fi3 

e33  _ - 


Cn,3333 


.(0) 


trial 


=  cne 


n+1  ~  u  n+ 1 

n(0)  -  c  f{0) 

A>*n+1  cnU(T /n+1 


ncn+ 1 


(3.32) 

(3.33) 

(3.34) 


(0) 

Cn+1  =  Cn 


(3.35) 


k  =  0 


(3.36) 


Step  2.  Update  stress  and  hardening  parameter: 

<ri+i  =  c„e!3,  - 

(3.37) 

«[%  =  a„  +  *A<*> 

(3.38) 

9n+l  —  '*■  U*n+1 ) 

(3.39) 

Step  3.  Check  for  failure  and  convergence: 

+  i 

II 

+  *; 

►-* 

+ 

+2 

I 

(3.40) 

If  ((lb  =  0  and  <  TOL0  or 

(  lvisc  =  false  and  |<?>l+i|  <  TOLi  and 

I17 l+i,33l  <  TOL2))  then  lquit  =  true 

If  (  lquit  =  true  )  then 

(k) 

<Tn+i  =  <rn+i 

(3.41) 

r  -  r{k) 

cn+l  —  cn+l 

(3.42) 

{k) 

ttn+l  —  Qn+1 

(3.43) 
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return 


endif 


Step  4.  Update  S A  and  e33: 


wS,  =  -  [d<r/!+i  •  nS,  +  w"  («!30; 


If  (  lvisc  =  false  )  then 

r  a<*+i>  i  r  <5a<*>  \ 

«s?  r 

wl+t 


,(*+n  = 

e33  J 


1  -1 


iVn+l,33 


33 


Cn,3333 


else 


-Ik)  T 

<n+1  =  “  '<*>  n^(‘! 


xX’iWi 

At(fc)  -  tk) 

4iIn+  1  —  ‘n+1 

(fc+1)  _  (fc)  .  (i) 
ln+l  ~  ‘n+1  +  ^‘n+l 


-At 

1  -  «*P  I  — 3T1 


n+1 


t 


n+1 


^  (<!»++i1)  -  <n+i  >  0)  then 


*tlnkh  =  a<£i  -  (d++i5)  -  *»+l) 


(k) 

<?n+ 1 
_(*) 
CTn+l,33 


lquit  =  true 


endif 


«'*+"  =  m'‘’  +  a<£,.v!2,/t 

(fc+1)  Cn,33iif]i  4-  Cn. 3322^22  +  2cni33i2ei2  +  2cn_2223^23  +  2 

£  - - 

Cn,3333 

rn.3333 
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(3.44) 

(3.45) 

(3-46) 

(3.47) 

(3.48) 

(3.49) 

(3.50) 

cn,3313^13 

(3.51) 


endif 


Step  5.  Update  return  direction  and  moduli: 


NL‘+,"  =  c  “UdafOli 


_(*+!)  _  -  f  gg  ^n+l 

cn+1  -  cn  0  A  ^ 


,‘~r‘  tvt(*+11 

cn+l  ‘  iNn+l 

Set  k  =  k  +  1  and  go  to  Step  2. 


(3.52) 

(3.53) 


Remarks: 

1.  Note  that  with  the  introduction  of  the  flag  “lvisc”,  Algorithm  7  can  be  employed  for 
both  elastic  and  rate- dependent  damage. 

2.  Because  6 A(fc+1)  is  explicitly  computed,  in  the  rate-dependent  case  —  see  (3.50)  —  the 

evaluation  of  =  e^/33  is  also  trivially  obtained  through  use  of  the  constitutive 

relation,  equation  (3.51). 

3.3  Applications  to  Plates  and  Shells 

The  above  algorithm  can  be  naturally  employed  in  the  analysis  of  structural  elements 
like  plates  and  shells.  An  excellent  source  of  information  concerning  three-dimensional 
continuum  based  shell/plate  formulations  is  Stanley[21].  Other  works  pertinent  to  the 
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present  subject  are  [16-18].  Because  of  the  completeness  of  [21],  only  major  features  of  the 
theory  relevant  to  the  present  discussion  will  be  recalled. 


1.  The  shell  element  is  obtained  by  degenerating  a  three-dimensional  element  with  lin¬ 
ear  displacement  variation  through  the  thickenss.  Two  assumptions  are  made  that 
distinguish  continuum  based  shell  elements  from  regular  three-dimensional  elements: 

Straight  Normals.  This  assumption  is  expressed  by  the  equation 

=  x(£,J7)  +  zx(f,7?)  (3.54) 

where  £  and  7/  stand  for  parent  domain  curvilinear  coordinates  while  z  is  a 
(linear)  through-thickness  coordinate.  We  note  that  z  is  usually  written  as 

*(6?)=(C-C)V2  (3.55) 

where  C,  €  [—1,1]  and  Q  locates  the  reference  surface,  h  is  the  shell  thickness. 
Inextensible  Normals.  This  can  be  described  incrementally  by  the  relations 


Au(£,77,;)  =  Au(£,t?)  +  zAu{^rj) 

(3.56) 

Au  •  x  =  0 

(3.57) 

However,  in  practice  it  is  necessary  to  maintain  the  inextensibility  condition 
exactly.  This  can  be  achieved  by  various  means.  For  example,  one  can  use  Au 
to  define  an  increment  of  angle  of  rotation  and  thereby  construct  an  orthogonal 
rotation  matrix  to  transform  the  normal  rigidly [21] .  or  one  can  use  a  radial 
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projection  algorithm  (e.g.,  Hughes  and  Liu  [7]). 


2.  An  orthogonal  local  coordinate  system  ^e^eLe^  referred  to  as  the  laminar  system 
is  introduced.  Its  most  important  attribute  is  that  e!3  is  normal  to  the  surface  (  = 
constant  at  each  point  in  the  current  configuration.  By  expressing  the  stress  tensor 
in  this  laminar  system,  the  constitutive  algorithms  introduced  above  for  the  plane 
stress  problem  can  be  easily  implemented. 

3.  The  continuum  based  approach  has  the  distinct  advantage  that  rigorous  finite  defor- 
mational  constitutive  theories  can  be  readily  utilized:  one  is  only  left  with  the  ZNS 
constraint  to  deal  with  as  described  above. 


5S 


Chapter  4 


Extension  to  Multiple  Failure 
Surfaces 

4.1  Rate- independent  Anisotropic  Elastic  Damage  Model 

Because  of  the  importance  of  models  that  need  to  accomodate  multiple  failure  surfaces 
(e.g.,  the  cap  model  [11]  or  those  proposed  by  Hashin  [4]),  we  present  extensions  of  the 
preceding  ideas  to  multiple  surfaces  [20,18]. 

A  theory  for  a  model  whose  failure  surface  is  defined  by  m  smooth  surfaces  is 
described  in  Box  4.  The  free  index  “A”  is  understood  to  take  on  the  values  1.2, ...  ,m. 
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Box  4.  A  multiple  failure  surface,  small-deformation  anisotropic  elastic  damage  constitu- 


(4.1) 

(4.2) 

(4.3) 

(4.4) 

(4.5) 

(4.6) 


Remarks: 


1.  H.  is  a  function  of  m  independent  strain-like  variables  ou,  a2 . om.  Thus,  the  total 

free  energy  may  be  written  as 

rp(e,c,aua2, . . .  ,am)  =  ^e  •  ce  +  H  (e*i, a2,  •  - . , am)  (4.7) 

2.  The  qA s  are  defined  by 

qA  =  -daAH  (4-S) 

3.  The  principle  of  maximum  dissipation  is  invoked  as  for  the  theory  of  Box  2.  In  the 
present  context,  this  leads  to 

ce  =  -  Y,  Aflc derfB  (4-9) 

s= 1 

=  \A  (4.10) 

4.  The  consistency  parameters  A'4  can  be  obtained  by  time  differentiating  the  damage 
functions  <t>A  and  setting  the  resulting  expression  to  zero.  From  (4.1),  (4.2),  (4.4) 
and  (4.10)  we  get 

4>a  =  d<rfA  ■  (ce  +  ce)  -  ]T  XBdlAaB'H  =  0  (4.11) 

B= 1 

Substituting  (4.9)  into  (4.11)  leads  to  the  final  system  of  equations  for  A'4, 

{A4}  =  [da fA  ■  cda/B  +  ^aeW]~l  {d*fB  •  ce}  (4.12) 


61 


5.  From  symmetry  arguments,  and  by  assuming  that  each  surface  contributes  a  rank- 


one  update  to  the  rate  of  change  of  the  elastic  moduli,  one  obtains 


B—\ 


cdcrfe  Q  cdg-  fB 
dcrfB  •  ce 


(4.13) 


4.1.1  Algorithm  for  the  Rate- independent  Anisotropic  Elastic 
Damage  Model 


We  consider  an  integration  scheme  similar  to  that  of  Algorithm  3;  however,  because  the 
active  surfaces  (i.e.,  those  for  which  <t> a  >  0)  are  not  known  a  prion ,  a  new  procedure 
needs  to  be  developed.  To  this  end  we  introduce  the  set  of  active  surfaces  at  iteration  k, 

4ct  =  {A  I  *!?’  >  0}  (4.14) 

The  procedure  is  then  defined  by  the  following  steps: 

1.  Let  ct  be  the  set  of  active  surfaces  at  the  fc-th  iteration.  Compute  increments 
6\A(-k\  A  €  ^act-  ky  employing  an  approach  similar  to  that  of  Algorithm  3.  as 
described  below. 

2.  Update  S\A^k^  by  setting 

6\A{k+1)  =  SXA[k)  +  (4.15) 

and  check  the  sign  of  <5A4(fc+!*.  If  negative,  drop  constraint  .4  from  the  active  set  and 
restart  the  iteration.  Otherwise  proceed. 
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The  algorithm  can  be  stated  as  follows: 


Algorithm  8.  Updated  explicit/implicit  algorithm  for  the  multiple  failure  surface,  small- 
deformation  anisotropic  elastic  damage  constitutive  theory  of  Box  J. 


Step  1.  Initialize: 


a^(0)  =  0 


_(o)  _  ^trial  _ 

&n+ 1  —  an+ 1  —  cnen+l 


i\d°)  _  c 

^/t,n+l  —  Cn'-,(T  J  A,n+ 1 


(0) 


c(0)  -c 

cn+l  —  cn 


k  =  0 


Step  2.  Update  stress  and  hardening  parameters: 


,n+l 


d”,  =  ^  -  E 
4.U.  =  o*.  +  a  s  4>t 

(k)  a  j{k) 

a4,n+l  ~  A 

(fc)  o  (  (*:)  (k)  Ik)  \ 

$A,n  + 1  —  (®l,n+l '  a2,n  +  l  i  •  ‘  ‘  >  ^m.n+l  J 


Step  3.  Check  for  failure  and  convergence: 


If  (k  =  0)  then 


C3 


(4.16) 

(4.17) 

(4.18) 

(4.19) 

(4.20) 


(4.21) 

(4.22) 

(4.23) 

(4.24) 


(4.25) 


jo)  _  AO)  (0) 

®A,n+ 1  ~  JA,n+l  +  9+.71+1  a A 


Jact  =  {'4  I  ®!Sn+i  >  °'  ^  = 


( 4. 20) 


<Z>+.n+l  =  /j^n+l  +  ^n  +  l  ~  <7^’  *4  ^  J'3 


endif 


If  ((Jfc  =  0  and  4Jt  =  0)  or  (|<^fc)n+1|  <  TOL  for  all  A  €  J^.))  then 


(fc) 

^  n+1  —  CTn+l 


(*) 

*-n+l  —  £-ti+1 


Ik)  4  _  r( 

aA,n+l  —  ®yl,n+l’  -i4  ^  a 


£*>1,71  +  1  -  Ofyl^,  >4  ^  J 


(4.28) 


(4.29) 


(4.30) 


(4.31) 


return 


endif 


Step  4.  Compute  £A-increments: 


For  .4,  B  eJi 


rtAB  (*)  _  / n  r(k)  -vH*)  i  -t+W  \ 

■^n+l  -  -  /.4,n+l  •  ^fl.n+l  +  doAcB  rtn  +  l  ) 


{6X^}  =  -  [z>£?1w]  {tfi+i} 

^temp  =  {-4  €  4't  I  >  0} 

^  (^act  ^  'Aemp)  ^hen 


(4.32) 


(4.33) 


(4.34) 


64 


(4.35) 

(4.36) 

(4.37) 


(4.38) 

(4.39) 


Remarks 


1.  We  note  that  <£^n+1  >  0  does  not  guarantee  that  surface  -4  will  ultimately  be  active. 
By  construction,  J^1)  C  J^ct* 

2.  The  linear  system  of  equations  (4.33)  reflects  a  Newton-like  iteration  to  force 

d>A,n+\  =0.  -4  €  Jact  (4.40) 
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The  similarity  between  (4.33)  and  (2.113)  should  be  noted. 


4.2  Rate-dependent  Damage  Model 


The  viscous  regularization  of  the  multi-surface  theory  follows  closely  the  development  of 
the  single  surface  theory  —  see  Section  2.5.  It  is  given  in  Box  5  for  completeness. 


Box  5.  A  multiple  failure  surface,  small- deformation  anisotropic  rate- dependent  damage 
constitutive  theory. 


Constitutive  equation: 


a  =  ce 


(4.41) 


Hardening  law: 


a,  —  _  Y'  (x  (d>g))  q2 

9-4 


(4.42) 


Damage  evolution  law: 


*  =  ■  E 


(X  (da))  c dtrfB  0  cd<rfB 
r  dafB  •  ce 


(4.43) 


6G 


4.2.1  Algorithm  for  the  Rate-dependent  Damage  Model 


Numerical  integration  of  the  equations  can  be  performed  the  same  way  as  for  the  single 
surface  case.  The  procedure  is  described  in  Algorithm  9  (cf.  Algorithm  5). 


Algorithm  9.  Updated,  explicit/implicit  algorithm  for  the  multiple  failure  surface,  small- 
deformation  anisotropic  elastic  damage  constitutive  theory  of  Box  4 ,  0T  rate-dependent 
damage  constitutive  theory  of  Box  5. 


Step  1.  Initialize: 


Ivisc  = 


true,  rate-dependent  damage  (Box  5) 
false,  elastic  damage  (Box  4) 


lquit  =  false 


A^n+l  —  tn+  j  —  tn 

(4.44) 

tnh  =  ^ 

(4.45) 

6Xa{0)  =  0 

(4.46) 

_(°)  _  -trial  _  _  , 

~n+l  °n+ 1  Cn^n-fl 

(4.47) 

T\J<°)  _  r  f)  f(°) 

+  l  ~  cnC'crjA,n  + 1 

(4.48) 

(0) 

Cn+1  =  Cn 

(4.49) 

k  =  0 

(4.50) 
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Step  2.  Update  stress  and  hardening  parameters: 


=  <4??*  -  E 


Be/ 


(*> 

act 


,<*) 

M,n 

- 


r(*) 


a>4,n+l  —  A  £  ^act 

fl(*)  __a  Wa(*>  ft<*>  Q(fc)  ^ 

Step  3.  Check  for  failure  and  convergence: 

If  (£  =  0)  then 


= fz» + tfu.  -  ^ 

Jf£t  =  {A  I  *2U.  >  0,  A  =  l,2,...,m} 


else 


<*<*>  _  fW  +  0<*>  4  f  iW 

YM.n+l  —  -M,n+1  ^  V^.n+l  ^  •'j 


act 


endif 


(4.51) 

(4.52) 

(4.53) 

(4.54) 


(4.55) 

(4.56) 


(4.57) 


If  ((fc  =  0  and  =  ®)  or 

(Ivisc  =  false  and  l^.n+il  5:  TOL  for  all  .4  €  ) )  then  lquit  =  true 


If  (  lquit  =  true  )  then 


( k ) 

^n+l  =  O-n+t 


Cn+1  =  C 


<*) 

n+1 


(fc)  1  7' 

a.4.n  +  l  —  <U„+I,  -4  €.  /&c^ 

6S 


r(fc> 


(4.5S) 

(4.59) 

(4.60) 


J  i d  7^9 

Q. 4,n  +  l  =  a4.n,  -*1  ? 


(4.61) 


return 


endif 


Step  4.  Compute  <5A-increments: 


For  A,B  £  4‘»t 


D^W  =  -  (fWj£+,  •  <1 +1  +  C.XA) 


If  (  lvisc  =  false  )  then 
{(5A',U+1)}  =  {<5A^<*>} 


D 


AB  W 

n+1 


1-1 


else 


xW  _  r 

*n+l  -  5 


a/W  _ 

^rn+l  —  ^n+ 1 


1  -  exp 


-At 


n+ 1 


;<fc) 


^71+  1 

(*)  ,  (*+l) 


=  a-1'1'  +  AtS.xiSl’./r,  a  e  j£'t 


(4.62) 


KU.}.  ABsC  (4.63) 


(4.64) 

(4.65) 

(4.66) 


endif 


Aemp  =  (-4  6  4ct  I  >  0} 


If  (^act  7^  1  tflen 


j(k )  _  r 
Jact  “  Jtem: 


For  .4  ^  .7 


(k) 

act 


6A-AfM  =  0 


(4.67) 


(4.6S) 


(4.69) 
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Go  to  Step  2. 


else 


If  (  lvisc  =  true  )  then 
(*+D  _  (JO  A  .(*) 

cn+l  ~  In+1  T  ^rn+l 

If  (C/1  -  Wi  >  0)  then 

A *L+i  =  AiS,  -  -  Wi) 

lquit  =  true 
endif 
endif 

t(*+1)  _  T 

Jact  _  "'temp 
endif 


Step  5.  Update  return  directions  and  moduli: 


N(*+i)  (*)  *  ,(*)  Ac  J 

“  Cn+ly(TJ/(,n+li  A  tz  J 


act 


.(*+D  _ 


-n+l 


=  cn  -  y.  6X 


) 

+  1 


act 


e  ■  N*fc+1) 

£n+l  i>/l,n  +  l 


Set  k  =  k  +  1  and  go  to  Step  2. 


(4.70) 


(4.71) 


(4.72) 


(.4.73) 

(4.74) 
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Remarks: 


1.  The  parameter  S  in  (4.64)  has  the  same  physical  meaning  as  —  v  'IX  \DoXh  in  the 
single  surface  algorithm:  cf.,  e.g.,  (3.46).  Because  we  now  have  a  matrix  DAB  to  deal 
with,  we  need  an  economical  way  to  bound  its  eigenvalues.  This  can  be  accomplished 
through  Gerschgorin’s  circle  theorem[15j.  We  define 


DaB  —  V(k)  nAB(k)  A  B  G  J 

u  —  X^.n  +  l-^n+l  >  A i  •D  t  J 


(fc) 

act 


(4.75) 


We  define  the  radii  of  the  Gerschgorin  circles  by 


Then, 


RA  = 


E  0 

B  €  4ct 
B  ^  A 


AB  | 


5  = 


min  [DAA  —  RA^j 


A  €  Jikl 


act 


(4.76) 


(4.77) 


This  value  of  5  provides  a  lower  bound  for  and  thus  provides  a  conservative 

estimate  for  stability  purposes. 


4.3  Plane  Stress  Algorithm  for  the  Rate-independent 
Damage  Model 


Generalization  of  the  single  surface  plane  stress  inviscid  algorithm  to  multiple  surfaces  is 
straightforward:  just  as  in  the  single  surface  algorithm,  we  introduce  an  extra  unknown 
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t33  and  an  extra  equation. 


(733  —  0  ( 4 .  /  b ) 

Derivation  of  the  algorithm  follows  along  the  same  lines  as  for  the  single  surface  algorithm 
described  in  Section  3.1.  It  is  presented  here  in  detail  for  completeness. 


Algorithm  10.  Updated  explicit/implicit  plane  stress  algorithm  for  the  multiple  failure 
surface,  small- deformation  anisotropic  elastic  damage  constitutive  theory  of  Box  4 ■ 

Step  1.  Initialize: 


a4(0)  =  o 

(4.79) 

€  =  e„+i 

(4.80) 

JO)  _  Cn.3311eU  +  ,3322 ^22  +  2c„i33i2€i2  +  2cn3323e23  +  2cn-3313e13 

(4.81) 

^33  — 

C«,3333 

_(°)  _  ,-trial  _  _  , 

^n+1  ^n+l  —  cn^n+l 

(4.82) 

T\j(°)  _  c  ft  /(°) 

^M.n+l  —  cnO(T  J  A  ,n+l 

(4.S3) 

r(°)  _  - 
tn+l  — 

(4.S4) 

k  =  0 

(4.85) 

Step  2.  Update  stress  and  hardening  parameters: 

=  c„c£,  -  £ 


(4. SG) 


oS+1  =  .-1  €  4ct  (4-S7) 

»St,  =  assi 

«!£U.  =  -a»,WKi+.."S+„---.aL‘,„+,)  (4.89) 

Step  3.  Check  for  failure  and  convergence: 

If  (k  =  0)  then 

4°)  _  r(°)  ,  (°)  _  /a 

<P.4,n+l  -  JA,n+ 1  +  <lA, n+1  -  aA  (4.9UJ 

Jact  =  {-4  I  >  0,  -4  =  1,  2, . . .  ,m|  (4.91) 

else 

tf.U  =  /S+i  +  «?,!.+.  -  <M  (4.92) 

endif 

If  ((k  =  0  and  j£>t  =  0)  or 

(kSjJil  <  TOLj  for  all  .4  €  j£t  and  l^i^l  <  TOL2))  then 
<rn+1  =  <t2,  (4.93) 

Cn+1  =  ch  + 1  (4.94) 

a^.n+i  =  a  If  L+ 1,4.  €  j£t  (4.95) 

Q4.n-H  =  Q.4.n,4  g  J&ct  (4.96) 

return 

endif 
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Step  4.  Compute  <5A-increments  and  €33  increment: 


For 

D™(k)  =  -  (dvfltLi  ■  N‘b‘L+1  + 

{6A'"‘+1i}  If  1 

ri  ^,.33  r 

.  ~~  { B.rx+ 1 .33/  cn.3333 

Jtemp  =  {A£  4ct  I  6xMt*"  >  °} 

IF  (-Fact  ^  ■Ftemp)  then 

r(^)  _  T 

•'act  —  •'temp 
For  .4  $  4‘>, 

6\A(k)  =  0 


(4.97) 


A,BiJ <‘>t  (4.98) 
(4.99) 


(4.100) 


(4.101) 


Go  to  Step  2. 

else 

^act  *  =  ^temp  (4.102) 

endif 
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Step  5.  Update  return  direction  and  moduli: 


J^rffc+ll  Ik)  p,  Ak)  ,  j(k) 

iN4.n  +  l  -  Cn+\U(rJA.n  +  n  A  =  Jact 


Ak+l 


“71  + 


1  =  Cn  ~  £  6X 


mm  iN^i1, 0  n^i', 


Aej 


(*+ 1) 
act 


.(*+i) 


T(*+l) 


en+l  '  ^/l,n+l 


Set  k  =  k  +  1  and  go  to  Step  2. 


(4.103) 

(4.104) 


4.4  Plane  Stress  Algorithm  for  the  Rate-dependent 
Damage  Model 


The  integration  algorithm  for  the  plane  stress,  multiple  surface  theory  with  regularization 
follows,  once  again,  along  the  lines  of  the  single  surface  theory.  At  the  end  of  the  inviscid 
calculation  the  viscoplastic  regularization  is  done,  the  33-strain  is  updated  to  account  for 
the  ZNS  constraint,  and  the  final  stress  is  computed. 


Algorithm  11.  Updated  explicit /implicit  plane  stress  algorithm  for  the  multiple  failure 


surface,  small-deformation  anisotropic  elastic  damage  constitutive  theory  of  Box  4,  or  rate- 
dependent  damage  constitutive  theory  of  Box  5. 


Step  1.  Initialize: 


lvisc  = 


true,  rate-dependent  damage  (Box  5) 
false,  elastic  damage  (Box  4) 


lquit  =  false 


&tn+l  —  t«+l  tn 

(4.105) 

t(°]  ~t 

ln+l  ~  cn 

(4.106) 

6A^(0)  =  0 

(4.107) 

c  =  cn+l 

(4.108) 

JO)  __  Cn.3311«ll  +  cn, 3322e22  +  2c„,3312ei2  +  2cn,3323e23  +  2cni3313ej3 

(4.109) 

^•33  — 

Cn,3333 

_(°)  _  -trial  _  _  , 

&n+ 1  —  **  n+ 1  —  cn^n+l 

(4.110) 

iv(0>  —  c  d*  f(0) 

iN4,n+l  —  cnU(TJA,n+ 1 

(4.111) 

JO)  _ 

Cn+1  —  Cn 

(4.112) 

O 

II 

-V 

(4.113) 
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Step  2.  Update  stress  and  hardening  parameters: 


rc+1  —  Cnen+1  /  *'iB,n  + 1 


BeJ 


(fc) 

act 


r(fc) 


^,U,=a4.»  +  «A'"‘>,.4€7<ct 

oji+i  =  “-4,.. -4  £  4ct 


Step  3.  Check  for  failure  and  convergence: 


If  (fc  =  0)  then 

^a.L+i  =  /a°i+i  +  IaI+i  ~  aA 

4ct  =  {4  I  ACh.  >  0.  4  =  1,2,.  ...m} 


else 


A*)  _  f(*)  ,  _(*) 

'Pa.n  +  l  —  JA,n+ 1  +  9a,n+l 


<7,4 ,  A  E  J. 


(*) 

act 


endif 


If  ((fc  =  0  and  =  0)  or 
(  lvisc  =  false  and  |<7>^„+i  I  <  TOL^  for  all  .4  G 
and  33 [  <  TOL2))  then  Iquit  =  true 

If  (  Iquit  =  true  )  then 

(k) 

&  n+  1  —  ^n+l 


(4.114/ 

(4.115) 

(4.116) 

(4.117) 


(4.1 IS) 

(4.119) 

(4.120) 


(4.121) 

(4.122) 


^  _  Jk)  ,  j(k) 

£*,4. n  +  1  ^.4.n  +  l’~*^  "'act 


(4.123) 


&.4. n  +  1  — 


(4.124) 


return 


endif 


Step  4.  Compute  <5A-increments  and  c33  increment: 


For  A,BeJ™ 


j-)AB  (fc)  _  (ft  t(k)  T>j(fc)  ,  fti  1 

Vn+ 1  -  -  \°cr  J a, n+1  •  ^B.n+l  +  aaA*B  rtn  + 1 ) 


(4.125) 


If  (  lvisc  =  false  )  then 


{a4<*+1>}  If  {6\A^}  1 


Jfc+D 

en+l,33 


Jk) 

en+l,33 


[■Oifi***]  t  ('«%) 

cn.3333 


°n+l,33  J 


e^rt  (4-J26) 


Ak)  _  l_ 
n+1  “  5 

A*<fe)  _?<*> 


A/1  '  —  p**  i  _  pvn 

1  ^n+1  ■*■  CXp 


MW+"  =  +  A («*,*!,*>  ,/r,  -4  €  Jit* 


— Atn+i 


(4.127) 


(4.12S) 


(4.129) 


cn,33Uell  +  cn. 3322^22  +  2cn,3312e13  4-  2cn3333f33  -f  2cn-33i3e13 


cn+l,33  ~ 


Cn, 3333 


cn,3333 


(4.130) 


endif 


^temp  ~  {-4  *=  ^act  I  ^  'M  !)  >  o| 


(4.131 


7S 


If  (Jact  7*“  ^temp)  then 


/(*•')  _  i 

Jact  ~  Jtemp 


(4.132; 


For  A  #  J[ 


6\A{k)  =  0 


(4.133) 


Go  to  Step  2. 


If  (  lvisc  =  true  )  then 


(fc+i)  _  ( k )  A  Ak) 
rn+l  —  rn+l  "*■  z-iIn+l 


(4.134) 


If  (<n+i  -  *n+i  >  0)  then 


a  =  a  -  (<SJi1)  -  Wi) 


(4.135) 


lquit  =  true 


endif 


endif 


t(*+1)  _  T 

Jact  —  •'temp 


(4.136) 


endif 


Step  5.  Update  return  direction  and  moduli: 


N(*+l)  (k)  ~  Ak)  j  r(K| 

+  l  *“n+l  J A.n  +  \ '  ^  *  3.Ct 


w"  =  C„  -  21  « 


g  ^A.n+l 
Jfc+«)  ivjffc+l) 

en+l  '  ^✓i.n  +  l 


( 4.13S ) 
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Set  k  -  k  +  1  and  go  to  Step  2. 


Chapter  5 


Conclusions 


In  this  report  we  have  considered  a  variety  of  three-dimensional  and  plane  stress  constitu¬ 
tive  models  and  algorithms  for  reinforced  concrete  plate  and  shell  structures.  Anisotropic 
damage  mechanisms  have  been  accounted  for  to  provide  a  setting  for  incorporating  vari¬ 
ous  failure  phenomena  within  a  homogenized,  or  distributed,  constitutive  representation. 
Rate-dependent  effects  have  been  introduced  by  way  of  a  viscous  regularization  technique. 
This  feature  is  useful  for  faithfully  modeling  high  rates  of  loading,  and  also  provides  a  con¬ 
stitutive  framework  which  avoids  certain  numerical  pitfalls  associated  with  strain-softening 
behavior.  Multiple  failure  surface  theories  have  also  been  investigated.  These  are  useful 
for  the  development  of  damage  theories  based  upon  often-used,  multiple  surface  theories, 
such  as  the  cap  model,  and  related  potentially  useful  theories  employed  in  the  modeling  of 
fiber-reinforced  composites. 
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